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Abstract 

The goal of this paper is to improve the performance of an electric motor by 
modifying the geometry of a specific part of the iron core of its rotor. To be more 
precise, the objective is to smooth the rotation pattern of the rotor. A shape 
optimization problem is formulated by introducing a tracking-type cost functional 
to match a desired rotation pattern. The magnetic field generated by permanent 
magnets is modeled by a nonlinear partial differential equation of magnetostatics. 

The shape sensitivity analysis is rigorously performed for the nonlinear problem by 
means of a new shape-Lagrangian formulation adapted to nonlinear problems. 

Keywords: electric motor, shape optimization, magnetostatics, nonlinear partial differ¬ 
ential equations. 


1 Introduction 

Advanced shape optimization techniques have become a key tool for the design of indus¬ 
trial structures. In the automotive and aeronautic industries, for instance, the reduction 
of the drag or of the noise are important features which can be reduced by changing the 
design of the vehicles. In general, when considering a complex mechanical assemblage, it 
is often possible to optimize the geometry of certain pieces to improve the overall perfor¬ 
mance of the object. In the industrial sector, the shape optimization of electrical machines 
is the most economical approach to improve their efficiency and performance. Shape op¬ 
timization problems are usually formulated as the minimization of a given cost function, 
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typical examples being the weight or the compliance for elastic systems. The most inter¬ 
esting and challenging problems of these type have linear or nonlinear partial differential 
equations constraints; see, for instance, [U 0 EH HU ESI HE1HH11221 EE1 123 EB1 ESI EH EE1 HD] 
and the references therein. 

In this work, a shape optimization approach is used to improve the design of an electric 
motor in order to match a desired smoother rotation pattern. As a model problem, we 
consider an interior permanent magnet (IPM) brushless electric motor consisting of a rotor 
(inner part) and a stator (outer part) separated by a thin air gap and containing both 
an iron core; see Fig. [I] for a description of the geometry. The rotor contains permanent 
magnets. The coil areas are located on the inner part of the stator. In general, inducing 
current in the coils initiates a movement of the rotor due to the interaction between the 
magnetic fields generated by the electric current and by the magnets. In our application, 
we are only interested in the magnetic field B for a fixed rotor position without any current 
induced. Since the magnetic properties of copper in the coils and of air are similar, we 
model these regions as air regions, too. We refer the reader to 0 EDI ED E2] for other 
approaches to the design optimization of IPM electric motors and to [43] for a special 
method of modeling IPM motors using radial basis functions. 

Due to practical restrictions, only some specific parts of the geometry can be modified. 
In our application, we identify a design subregion fl of the iron core of the rotor subject to 
the shape optimization process. Our objective is to modify fl in order to match a desired 
rotation pattern as well as possible. Practically this is achieved by tracking a certain 
desired profile of the magnetic flux density, which is done by reformulating the problem 
as a shape optimization one by introducing a tracking-type cost function. 

The shape optimization of fl has been considered in m from the point of view of the 
topological sensitivity [34, 39j. However, the derivation of the so-called topological deriva¬ 
tive for nonlinear problems is formal since the mathematical theory for these problems is 
still in its early stages; see HHQ1 ED for a few results in this direction. Moreover, the 
drawback of the topological derivative is that it usually creates geometries with jagged 
contours. 

In this paper, we focus on the shape optimization of the design domain f2 by means of 
the shape derivative , which, contrarily to the topological derivative, proceeds by smooth 
deformation of the boundary of a reference design. In this way, the optimal shape has a 
smooth boundary provided that the numerical algorithm is carefully devised. Computing 
the shape derivative for problems depending on linear partial differential equations is a 
well-understood topic; see for instance [ZliniHD!- For nonlinear problems, the literature 
is scarcer and the computation of the shape derivative is often formal. A novel aspect of 
this paper is to provide an efficient and rigorous way to compute the shape derivative of 
the cost functional without the need to compute the material derivative of the solution 
of the nonlinear state equation. The method is based on a novel Lagrangian method for 
nonlinear problems and on the volume expression of the shape derivative; see [23 SI]. 
This allows to obtain a smooth deformation field used as a descent direction in a gradient 
method. In the numerical algorithm, the mesh is deformed iteratively using this vector 
field until it reaches an equilibrium state. 

The rest of the paper is organized as follows: In Section [2j we formulate the shape op¬ 
timization problem and give the underlying nonlinear magnetostatic equation. Existence 
of a solution to the shape optimization problem is shown in Section [3| In Section [4| we 
introduce the general notion of a shape derivative and give an abstract differentiability re- 
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(a) (b) 

Figure 1: (a) Geometry of the electric motor with permanent magnet areas Q ma (or¬ 
ange with arrows indicating the directions of magnetization), ferromagnetic material Q T f 
(brown), coil areas fl c (light blue), air Q a (dark blue) and thin air gap (yellow layer be¬ 
tween rotor and stator); r 0 is a circle located in the air gap. (b) Zoom on the upper 
right quarter of the motor where the reference design area f2 ref C corresponds to the 
highlighted (red) region. The reference region Q ref is subject to the shape optimization 
procedure. 


suit which is used later on to compute the shape derivative of the cost functional. Section 
[5] deals with the shape derivative of the cost functional. Finally, in Section [6| a numerical 
algorithm is presented to optimize the design of Q, and numerical results showing the 
optimal shape are presented. 


2 Problem formulation 

Let D C R 2 be the smooth bounded domain representing an IPM brushless electric motor 
as depicted in Fig. [l] with a ferromagnetic part Qj cf C D. permanent magnets fi rna C D, 
air regions fl a C D and coil areas Yl c . The design domain Q is included in a reference 
domain Q ref c Oy ef . The inner part of the motor is called the rotor and the outer part 
the stator. They are separated by a small air gap, the thin yellow circular layer in Fig. 
[l} By To we denote a circle within the air gap. When an electric current is induced in 
the coils, the rotor containing the permanent magnets rotates. In reality the motor is 
a three-dimensional object, but considering the problem only on the cross-section of the 
motor is a modeling assumption that is commonly used; see U E]- For a comparison 
between two- and three-dimensional models of electric motors, see [251 S2j- 

Denote T := dVt the boundary of the optimized part Yl which is assumed to be Lips- 
chitz. We introduce the variable ferromagnetic set Ylf := (0) uf \ O ref ) U Yl and Yf := dYlf. 
The permanent magnets create a magnetic field in D. In our application, we assume the 
coils to be switched off. Thus, no electric current is induced and the rotor is not moving. 
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The magnetic field generated by the permanent magnets can be calculated via a boundary 
value problem of the form 


— div(/?n(x, | Vu| 2 )Vu) = / in flf and D \ Qf, 
u = 0 on dD , 

with the transmission conditions on the interface T f 

M = 0 on r /> (2) 

\Pn{x,\Vu\ 2 )d n uj = 0 on T f , 

where n denotes the outward unit normal vector to Qf. Defining v~ the restriction of 
some function « on SI/ and v + its restriction on D \ Qf we denote by [u] the jump of v 
across the interface Tf, i.e. 

M = v+ \r f ~ v~\v r 

The nonlinear, piecewise smooth function f3 is defined for (x, () G D x R as 

Pn(x, C) := Pi(C)xn f (x) + P2(0x D \Uj( x ) 

= Pl(0(xn( x ) + Xn,ref\tt^f( X )) + P2(0(.X D \^ref( X ) + X^ef\n(^)) , 

where x is the indicator function of a given set. Note that the expression above is mean¬ 
ingful since Q C Q ref C The weak form of 0 reads 

Find u E Hq(D) such that / Pq(x, \ Vu| 2 )Vu • Vudx = (/, v) Vv E Hq(D), (3) 

Jd 

where (•, •} denotes the duality bracket between H~ 1 (D ) and Hq(D). The scalar function 
u is the third component of the vector potential of the magnetic flux density in three 
dimensions, B = curl((0, 0, u) T ). In our model, we consider the restriction of B to a 
two-dimensional cross-section since the third component vanishes. 

In the sequel, we make the following assumption for /3 1 and /? 2 : 

Assumption 1 . The functions (3i , /? 2 : R —» R satisfy the following conditions: 

1. There exist constants ci, 02 , 03,04 > 0, such that 

ci < Pi(0 < c 2 , c 3 < fi 2 ( C) < c 4 for all (GR. 

2. The function s i-» Pi(s 2 )s is strongly monotone with monotonicity constant m and 
Lipschitz continuous with Lipschitz constant L: 

(Pi(s 2 ) s — Pi(t 2 ) t ) (s — t) > m(s — t ) 2 for all s, t > 0, 

| Pi{s 2 ) s — Pi(t 2 ) £)| < L|s — t| for all s,t > 0 . 

3. The functions Pi,p 2 ore in C , 1 (R). 

4 . There exist constants A, A > 0 such that for i = 1,2, 

AM 2 < Pi(\p\ 2 M 2 + mp\ 2 )(v ■ P) 2 < A|h | 2 for all V , p G R 2 . 
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The task is to modify the shape of the design region Cl C fi ref in such a way that the 
radial component of the resulting magnetic flux density on the circle To in the air gap fits 
a given data as good as possible. We consider the following minimization problem: 

{ minimize J(Cl,u ) := / \B r {u) — Bd\ 2 ds 

dr 0 (4) 

subject to Cl G O and u solution of i§ 

where 

O = {fl C fd ref C Cljf Cl open and Lipschitz with uniform Lipschitz constant Lo} (5) 

and T 0 C D is a smooth one-dimensional subset of D. The sets Q rcf and Cl T f are reference 
domains; see Fig. [I] Here, B r (u) = Vu ■ r where r is the tangential vector to T 0 and 
B,[ G C' 1 (T 0 ) denotes the given desired radial component of the magnetic flux density along 
the air gap. In order to obtain the first-order optimality conditions for this minimization 
problem we compute the derivative of J with respect to the shape Cl. 

Remark 1 . Let us note that 


• in our application, the right-hand side f represents the magnetization of the perma¬ 
nent magnets. In general, it can he a combination of magnetization and impressed 
currents in the coils. 


• the link between fi and the magnetic reluctivity is fi(s) = u(y/s). In this case, the 
boundary value problem 0 is called the two-dimensional magnetostatic boundary 
value problem; see \Wf . We will see in Section [d| that Assumption^ 7] is satisfied 
for v due to physical properties. 

• the Dirichlet condition u\qd = 0 implies B ■ n = 0, thus no magnetic flux can leave 
the computational domain D. 


Given Assumption [lj we can state existence and uniqueness of a solution u to the 
state equation ([3j). 

Theorem 1 . Assume that Assumption 7|7 and 7|1 hold. Then problem Q admits a 
unique solution u G Hq(D) for any fixed right-hand side f G H~ 1 {D ) and we have the 
estimate 

1 \\ h 1 ( d ) < cH/Htf-ip). 


u 


Proof. A proof can be found in |T6, 35, 44j. 


□ 


3 Existence of optimal shapes 

In this section, we prove that problem s has a solution O*. We make use of the following 
result m Theorem 2.4.10] 

Theorem 2. Let Cl n be a sequence in O. Then there exists Ll G O and a subsequence 
Ll nk which converges to Cl in the sense of Hausdorff, and in the sense of characteristic 
functions. In addition, fl nk and dCl nk converge in the sense of Hausdorff respectively 
towards and dCl. 
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Let fl n G O be a minimizing sequence for problem Q . According to Theorem [2j we 
can extract a subsequence, which we still denote Q n , which converges to some 0 <S O. 
Denote u n and u the solutions of with 0 n and Q, respectively. We prove that 

u n —y u in H ( \(D). 

First of all in view of Theorem Q] we have 


11 11 #!(£>) < c\\f\\ H ~'( D ) ( 6 ) 

where c depends only on D. Thus we may extract a subsequence u n such that u n —>■ u* in 
L 2 (D ) and Vi/.„ —^ Vu* weakly in L 2 (0). Extracting yet another subsequence, we may as 
well assume that —y IT in the sense of characteristic functions applying Theorem [2j 

Subtracting the variational formulation for two elements u n and u m of the sequence 
and choosing the test function v = u n — u m we get 


/ ((/?n„(®, |Vu n | 2 )Vu n - (/3n m (x,\Vu m \ 2 )Vu m ) ■ V(u n - u m ) dx = 0 

Id 


(7) 


Let us introduce for simplicity the notation /?" := /?i(|Vu n | 2 ) and (3% := /3 2 (\ Vu n | 2 ). Then 
0 becomes 


' D 


(/?i (x« n + Xnf ! \nret) + @2 (Xr>\n r , ef + Xn-\nJ)V«n • V(w n — 


U r 


~ (PTiXfim + Xn^ ef \n ref ) + P 2 (XD\n r f + Xn ref \n m ))Vu m • V(w n — u m ) dx — 0 

This yields 


/] | /-J : /.{ ; /.( | /r, •• /(> — 0 


( 8 ) 


where 


7i := / (xn„ - XoJ/^Vun • V(« n - u m ) dx , 


' D 


h := / Xn m (/3rVu n - /5rVu m ) • V(u n - u m ) dx, 


ID 


h := / Xn-f\n-f(^Vu„ - • V(« n - u m ) dx, 


/£> 


I 4 . / XD\n r f(02 ^ u n $2 m ) ' U m ) dx. 


ID 


h := / (Xn«f\nn “ * V(u n - u m ) dx, 


ID 


Iq . / (/?2 ^U n @2 * V(' 2 / n Um) dx. 

JD 

To estimate the above integrals, we use the following lemma: 

Lemma 1. Letp,q E R 2 . If Assumption 1^2 holds, then 

m\p\ 2 )p - A(|?l 2 )?] • (p - 9) > m|p - 
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Proof. Using Assumption [Tp] we compute 

\PM 2 )p - PM 2 )q] ’(p-v) 

= m\p - q\ 2 + [(A(|p| 2 ) - m)p - (AO#) “ rn)q] • (p - q) 

= m\p - q\ 2 + (A(|p| 2 ) - m) ( \p\ 2 - p ■ q) - (AO#) - m) (p ■ q - \q \ 2 ) 

>IpI(Ip|-M) >° <M(IpI-M) 

= m\ P -q\ 2 + m P \ 2 )\p\ - m\ 2 )m\p\ - i# - i# 2 

S -v-' 

>m{\p\-\q \) 2 

> m\p — q | 2 . 


□ 


Applying Lemma [l] with p = Vu n and q = V« m we get 


I 2 + ^3 + h + I6 > m / (%n m + Xn ref \nref + XD\ci T f ei + Xn TeS \n m )\Vu n — \7u m \ 2 dx 


id 


= m \Vu n — Vw m | 2 dx. 


(9) 


ID 


Holder’s inequality yields 


\h\ < ||A l ||z,°°(D)||xa l — Xn m ||L i -(D)||Vn n ||i / 2( jD )||V(n n — n m )|| L qx)) 

with r _1 + s -1 = 1/2, r,s > 1. Performing a similar estimate for / 5 and in view of 
Assumption | 1 | 1 | this yields 

\h\ < c 2 ||xn n — XnJ\Lr { D)\\Vu n \\ L 2 {D) \\\7(u n — u m )|| £,*(£>) (10) 

\h\ < c 4 ||xn„ — Xn m |U’'(D)||V'U n ||i < 2 ( £ ))||V('U n — u m )\\ L s^ D ) (11) 

Using equality (| 8 ]) and inequalities (|6|),(|9l),([l0|),([IT|) we obtain the estimate 

/ |Vw n — Vtq# dx < c\\xn n — Xn r J|Lq#|/||rr- 1 (.D)||V('U n — u m )|| £,«(£>) 

Jd 

Since xci n is a characteristic function, the parameter r can be chosen arbitrarily large, 
and consequently, s can be chosen arbitrarily close to 2. Therefore, assuming a little more 
regularity than H 1 for the solution of 0 . the convergence of the characteristic functions 
of Cl n in L P (D) implies the strong convergence of u n towards u* in Hq(D). Consequently, 
we obtain the following result: 

Proposition 1. Let Ll n G O be a minimizing sequence for problem 0 and Ll be an 
accumulation point of this sequence as in TheoremAssume there exists e > 0 such that 
the solution u of Q satisfies 

IMIif 1 + £ (D)niT ( h.D) A c 

where c depends only on f and D. Then the sequence u n G Hq(D) corresponding to Cl n 
converges to u strongly in Hq(D), where u is the solution of 0 in Ll. 
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Proof. We have seen that there exists u* E H ( \(D) such that u n —y u* in H { \(D). We 
just need to prove that u* = u. The strong convergence of u n in H t \ (D) yields Vu„ —>■ 
V//* pointwise almost everywhere in D, and also the pointwise almost everywhere (a.e.) 
convergence (3™ —> /?,;(|Vi/*| 2 ) for i = 1,2. We also have the pointwise a.e. convergence 
of the characteristic function xci n 1° Xn* which implies the pointwise a.e. convergence 
/Jn n (-, |V u n | 2 ) — > /3q*(-, |Vu*| 2 ). Next, the weak formulation for u n is 

[ fin n (x, |Vu n | 2 )Vu n • Vvdx = (f,v) for all v E Hq(D), 

Jd 

The strong convergence of u n in H ( \(D) and the pointwise convergence of 3n„ (•, Vu n | 2 ) 
implies 

[ p n *(x, |Vn*| 2 )Vn* • Vvdx = (f,v) for all v E 

Jd 

which proves finally u* = u. □ 

Remark 2. In fact we have proven a stronger result, i.e., the Lipschitz continuity of u n 
in Hq(D) with respect to the characteristic function xn n - 


4 Shape derivative 

4.1 Preliminaries 


In this section, we recall some basic facts about the velocity method in shape optimization 
used to transform a reference shape; see mug. In the velocity method, also known as 
speed method, a domain fl C D C R 2 is deformed by the action of a velocity field V 
defined on D. Suppose that D is a Lipschitz domain and denote its boundary E := dD. 
The domain evolution is described by the solution of the dynamical system 


d 

dt 


(t) = V(x(t)),t E [0,e), 


x(0) = X eR 2 


( 12 ) 


for some real number e > 0. Suppose that V is continuously differentiable and has 
compact support in D , i.e. V E V l (D, R 2 ). Then the ordinary differential equation (12) 
has a unique solution on [0, e). This enables us to define the diffeomorphism 


T t : R 2 -> R 2 ; X i—^ T t (X) := x(t). (13) 

With this choice of V, the domain D is globally invariant by the transformation T t , i.e. 
T t (D ) = D and T t (dD) = dD. For t E [0, e), T t is invertible. Furthermore, for sufficiently 
small r > 0, the Jacobian determinant 


m := det DT t (14) 

of T t is strictly positive. In the sequel, we use the notation DT p 1 for the inverse of DT t 
and DTff T for the transpose of the inverse. We also denote by 

Ut) ■■= m\DT t ~ T n\ (15) 

the tangential Jacobian of T t on dD. 

Then the following lemma holds [7j: 
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Lemma 2. For ip G W^(R 2 ) and V" G P 1 (R 2 ,R 2 ) we have 

V(<y? o T t ) = DTj(Vp) o T t , o T t ) = (Vp • V) o T t , 

= £(i) [divH(t)] o T t , Cr(0) = div 9£ > V := divV\ aD - DVn ■ n. 

Definition 1. Suppose we are given a real valued shape function J defined on a subset 
H of the powerset 2 R . We say that J is Eulerian semi-differentiable at ft G S in the 
direction V if the following limit exists in R 

<17(11; 10 := Urn«“>. 

K J t\ 0 t 

If the map V —x dJ(i 2; V) is linear and continuous with respect to the topology of 
V(D, R 2 ) := Cf°(D, R 2 ), then J is said to be shape differentiable at fl and dJ{fl\V ) 
is called the shape derivative of J. 


4.2 An abstract differentiability result 

Let E and F be Banach spaces. Let G be a function 

G:[0,t]x£xF4R, (t, ip, xp) i—x G(t, p>, xp) (16) 

and define 

Eft) := {u e E\ dipG(t,u,0]xp) = 0 for all xp e F}. (17) 

Let us introduce the following hypotheses. 

Assumption 2 (HO). For every ( t,xp ) G [0, r] x F. we assume that 

(i) the set Eft ) is single-valued and we write Eft) = {if}, 

(ii) the function [0,1] 3 s i—x G(t, sxfi + sfa 1 — u°),xp) is absolutely continuous, 

(Hi) the function [0,1] 3 s e-X d v Gft,su l + (1 — s)u°,xp]ip) belongs to L 1 (0,1) for all 
<peE, 

(iv) the function xp h- x G(t,u,xp) is affine-linear. 

For t G [0,r] and xfi G Eft), let us introduce the set 

Y ft, xf, u°) := G F| \/p G E : J d v Gft, su f + (1 — s)u°, q ; <p) ds = oj- , (18) 

which is called solution set of the averaged adjoint equation with respect to t , xf and u°. 
Note that Y(0,u°,u°) coincides with the solution set of the usual adjoint state equation: 

Y (0, u°, u°) = [q G F\ d v G( 0, u°, q;ip) = 0 for all (p G E} . (19) 

The following result, proved in m allows us to compute the Eulerian semi-derivative of 
Definition [l] without computing the material derivative u. The key is the introduction of 
the set m. 
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Theorem 3. Let Assumption (HO) hold and the following conditions be satisfied. 

(HI) For all t G [0, r] and all ip G F the derivative dtG(t, u°, exists. 

(H2) For all t G [0, t], the set Y(t, vf, u°) is single-valued and we write Y(t, vf , u°) = {p f }. 

(H3) For any sequence of non-negative real numbers (t n ) ne n converging to zero, there 
exists a subsequence (£ nfc )fceN such that 


lim d t G(s, u°,p tn k) = d t G(0, u°,p°). 

k->-oo 

s \0 


Then for if G F we obtain 


j t (G{t,u\fi)) | t=0 = d t G(0,u°,p°) 


( 20 ) 


4.3 Adjoint equation 

Introduce the Lagrangian associated to the minimization problem © for all <p, ip G Hq (D): 

G(Ll, ip, if) := f \B r (ip) — B d \ 2 ds + f fi(x, | Vw| 2 )V(p • Vxp dx — (/, ip) (21) 
J r 0 Jd 

The adjoint state equation is obtained by differentiating G with respect to <p at Lp = u 
and if = p, 

d v G(Cl,u,p; ip) = 0 for all ip G Hq(D), 

or, equivalently, 

2 / d^/3(x, | Vu\ 2 )(Vu ■ V<p)(Vw ■ Vp) dx + / fi(x, |Vw| 2 )Vp • Vpdx 
Jd Jd (22) 

— —2 ( B r {u ) — Bfi)B r {<p) ds for all <p G Hq(D). 

J r 0 

Introduce the mean curvature k of T 0 and the Laplace-Beltrami operator A T on T 0 : 

du d 2 u 

A t u := A u - k-- — 7. 

an on z 

Using B r (u) = V r u as well as the equalities 

(Vu ■ V<p)(Vw ■ Vp) = ((Vu <8> V«)Vp) • Vp, 

/ (' V T u — Bfi) • V T p ds = — A T mp — pnBd ■ n ds 

Jt 0 Jt 0 


and Green’s formula, we deduce the corresponding strong form of (22) 


— div(Ai(Vw)Vp) = 0 in 

— div(.4.2(Vu)Vp) = 0 in D\Q, 

p = 0 on dD, 


( 23 ) 
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with the transmission conditions 


b] r = ° ° n r > 

[A{F7u)F7p ■ n\ r = 0 on T, (24) 

[Aiyu)\7p • n] ro = 2 (A T u — uB d ■ n) on r 0 , 

where 


A(Vu) := A(Vn)xn + A 2 (\7u)xd\q, 

Ai(Vu) := (3i(\ Vm| 2 )/ 2 + 2d^/3 i {\Vu\ 2 )'Vu 0 Vu £ R 2,2 , * = 1,2. 


Note that, with this notation, the variational form of the equation can be written as 


/ A(Vu)Vp ■ Vp dx =—2 / (B r I 
J d J r 0 


u 


B d )B r {p) ds for all p £ Hq(D). (25) 


Now let us investigate the existence of a solution for the adjoint equation 
Theorem 4. Let Assumption ^f^ hold. For given u £ Hq{D) the equation 

AfS7u)F7p ■ Vpdx = — 2 / (B r (u) — B d )B r {p) ds for all p £ Hq{D). 


(26) 


Jd Jr 0 

has a unique solution p £ Hq(D). 

Proof. For fixed u £ H t \ (D). define the bilinear form 


a\u- •) : Hq(D) x Hq(D) R 

(v,w) (->• / A(Vu)V v ■ V w dx. 

Jd 

We check the conditions of Lax-Milgram’s theorem. The ellipticity of the bilinear form 
a'{up,-) can be seen as follows: 

a'{w,v,v) = / A(Vu)Vv ■ Vv dx > \ / |Vu| 2 da; > A C'||v||#i( D ) 

Jd Jd 

where we have used the first estimate in Assumption EH and Poincare’s inequality since 
v £ Hq{D). The boundedness of the bilinear form a'(up,-) can be seen by Hblder’s 
inequality and again Assumption EB The right-hand side is obviously a linear and 
continuous functional on Hq(D), 


(F u , p) = -2 f ( B r (u ) - B d )B r {p ) ds, 

J r 0 

thus the theorem of Lax-Milgram yields the existence of a unique solution p to the varia¬ 
tional problem 


a!{u; p,p) = ( F u , p) for all p £ H^{D). 


□ 
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5 Shape derivative of the cost function 

In this section we prove that the cost function J given by is shape differentiable in the 
sense of Definition[l} Moreover, we derive a domain expression of the shape derivative. To 
be more precise, Theorem [3] is applied to show Theorem [5} Anticipating on the application 
of Section [6j we assume in what follows that / has the form 

/ = /o + divM with M = + M 2 XD\n ma ( x ) 

where /o G if 1 (D). 

In this section we assume Ll C O ref , V G XT(R 2 ,R 2 ) and supp(V) fi To = 0. Denote 
Lljff k = 1,8, the connected components of Q ref (see Fig. [lj. Introduce T ref the 
boundary of Q ref . The four sides of r ref are denoted r^ ef,JV , r^ ef ’ , T). ef ’ E , r^’ 5 where the 
exponents mean “north”, “south”, “east”, “west”, respectively. We assume V ■ n = 0 
on F k ei ’ s and V ■ n < 0 on r(. cf,i: U T^, ef,VF U F( rf:A? . These conditions guarantee that 
Q t := T t (Q) c D ref . In addition, we assume that the vector held V is such that the 
transformation T t satishes T)(Q rna ) = f2 ma for t small enough. 

Theorem 5. Let /A and (3 2 satisfy Assumption [7} Then the functional J is shape differ¬ 
entiable and its shape derivative in the direction V is given by 


dJ{fl- V) = — j (/ 0 div(V) + V/ 0 • V)pdx 


ID 


+ f P'(O) Vp • Mi dx + f 

J ^ma J X 

+ [ Pn( x , |Vu| 2 )Q'(0)Vm • Vpdx 


• Mi dx + / P ; (0) Vp • M 2 dx 

D\£} m a 


(27) 


id 


L 


— / 2d^/3a(x, | \7u\ 2 ){DV 1 Vw • Vu)(Vu • Vp) dx 


where P'(O) = (div V)I 2 - DV T , Q'(0) = (div V)I 2 - DV T - DV, I 2 G M 2 ’ 2 is the identity 
matrix, and u,p G Hq(D) are respectively the solutions of the problems ((Tj) -((2|) and (23)- 

Remark 3. Note that the last integral in (27) is well-defined thanks to Assumption [ip[ 
To see this, note that V G C 1 (D, R 2 ), and that, for all f G R 2 , we have / d / (lC| 2 )lC | 2 < A. 
Hence 


Id 


2d c (3 n (x, \Vu\ 2 )(DV T Vu ■ Vu)(Vu • Vp) dx 


<C d ( (3 n (x,\Vu\ 2 )\Vu\ 2 \Vu-Vp\dx 


id 

< CA / |Vu • Vp| dx < oc. 

Jd 


(28) 


The other terms in (|27j) are obviously well-defined. 

Proof of Theorem [5j Let us consider the transformation T t defined by (13) with 
V G T) l (D, R 2 ). In this case, T t (D ) = D, but, in general, T t (f2) := Ll t fl. We define 
the Lagrangian G(Ll t , (p,ip) at the transformed domain Q t for all <p, fi in H ( \ (D): 


G(fi*,<p,V0= / \B r (ip) - B d \ 2 ds + / Mx,\Vp\ 2 )Vp-V^da;-(/,^). 
JYo Jd 
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Since 


/ = f 0 + div M with M = + M 2Xd\dE~S x ) 

where M\ and M 2 are constant vectors, we transform the last term in G to 

(/,'*/’)= [ M + (div M, ip) = f fo'ip f M ■ VV>, 

Jd Jd Jd 

which yields 

G(fl t ,P,^) = f \B r ((p) — B d \ 2 ds + f /3n t (x, |V<p| 2 )V<p • da: 

7r 0 

/ /ot; I f A#! • V'0 + f _M 2 -V^. 

J D J O ma J D\£} ma 

In order to differentiate (/?, with respect to £, the integrals in (/?, *0) need to 

be transported back on the reference domain ^ using the transformation T t . However, 
composing by T t inside the integrals creates terms ip o T t and if; ° T t which might be 
non-differentiable. To avoid this problem, we need to parameterize the space Hq(D) by 
composing the elements of Hq(D) with T t _1 . Following this argument, we introduce 


® (t, tp) := V°T t 1 ,ipoT t 1 ). 


In (29), we proceed to the change of variable x = T t (x). This yields 


(29) 


®(t,ip,xp) = £r(t)\B r ((p) — B d o T t \ 2 ds 

JT 0 


+ [ Pci(x, |M(t)V<^| 2 )M(t)V</? • M(f)V^>£(t) dx 

Jd 

- j f 0 o T t ^(t) dx+ f M l - M(t)Vxp^t) dx 

J D J O m a 

+ [ _ M 2 -M(t)Vip£(t)dx, 

J D\ Dma 


(30) 


where M(f) = DTf r and £(t),£ r (£) are defined in (14) and (15), respectively. Note that 
we have used the assumption T t (Q ma ) = f2 ma in the computation of ®(f, p,ip). 

Note that J(f2 t ) = <&(t, u t ,xp) for all xp G Hq(D), where vJ G Hq(D) solves 


[ p n (x, |M(f) W| 2 )M(t) W • M (t)Vip£(t) dx 
Jd 

= f fo° T t xp£(t) dx- f Mi • M (t)Vxp£(t) dx - 

J D j rima 


'D\n n 


M 2 ■ M(i)vv>f (t) dx. 


(31) 


To prove Theorem [5j we need to check the conditions of Theorem [3] for the function 
<&(t, xp, ip) with E = F = Hq(D). 

Verification of (HO). Condition (H0)(i) is satisfied since E(t) = {«*}, where vJ G 
Hq(D) is the solution of the state equation ( |M| ). Conditions (ii) and (iii) of Assumption 
(HO) are also satisfied due to the differentiability of the functions /?i,/? 2 and Assumption 
1. Condition (H0)(iv) is satisfied by construction. 
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Verification of (HI). Condition (HI) is satisfied since M (t), £(t) and £ r (t) are 
smooth. 

Verification of (H2). Yit^u 1 ,u°) = {p*}, where p l G Hq(D) is the unique solution 
of 



Zdtfsiix, |M(t)V«?| 2 )((M(t)V«? <g> M(t)Vuf)M(t)Vp*) • M (t)Vip)£(t) dxds 


0 JD 



/0 JD 



Pn(x, |M(t)V<| 2 )M(t)V^ • M(t)Vp^(t) dxds 

2 (B r {u s t ) — Bd)B r (pld)^(t) dxds for all ^ G Hq(D). 


+ 


>o J r 0 

which can also be rewritten in a more compact way as 

«4(M(t)V<)M(t)Vp* • M(t)V^(t) dxds 

f f 2(B r (u s t ) — Bd)B r (ip)£(t) dx ds. 
o Jr 0 ' ' 


(32) 



0 J D 


(33) 


To prove that the previous equation has indeed a unique solution, we first check that all 
integrals are finite in the previous equation. To verify this, we use Holder’s inequality to 
obtain 


id 


2d c p n (x, |M(t)Vu?r)((M(t)V< <g> M(t)V<)M(t)Vp*) • M(t)V^)f(t) dx 


< C (^J 2 d c fa(x, |M(t)V<| 2 )(M(t)V< • M(£)Vp*) 2 dx'j ■ 
J 2d < f3 n (x, |M(t)V<| 2 )(M(t)V«? • M(t)V^) 2 dx'j 


and 


id 


Pn(x, |M(t)V<| 2 )M(t)V^ • M(t)Vp^(t) dx 


< c 


’D 


\ 1/2 

fci(x, |M(t)V<| 2 )|M(t)V^| 2 dx j 


\ 1/2 

fci(x, |M(t)V<| 2 )|M(t)Vp*| 2 dx 


id 




Adding both equations and using part [4] of Assumption [TJ we get 

_4(M(t)Vti))M(i)Vp‘ • M(t)V^(t) dxds < c||</’IIh‘(B)IIJ j, IIh‘(b). 


(34) 


id 


where the constant c > 0 is independent of s. 

The existence of a solution p t follows from Theorem [tJ Since M(i) = DT t T , there are 
numbers C > 0 and r > 0 such that, for all t G [0, r] and p G R 2 , we have M (t)p-p > Cjp| 2 . 
Note that p° = p G V(0, u°, u°) is the unique solution of the adjoint equation (|23|)-(24). 


Verification of (H3). To verify this assumption we show that there is a sequence 
(p tk )ke n, where {p tk } = Y{tk,u tk ,u°), which converges weakly in Hq(D) to the solution 
of the adjoint equation and that (t,ip) h-» 0 t &(t, u°, b) is weakly continuous. In order to 
prove this, we need the following lemmas. 
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Lemma 3. Let m E {0,1} and the velocity field V E T> 1 ( R 2 ,R 2 ) be given and p E 
H m ( R 2 ). We denote by T t , the transformation associated to V. Then we have 

lim p o T t = ip and lim p o Tfi 1 = p in H m (D). 


Proof. See for instance [7J. 

Recall that according to Remark [2] there are constants C,e > 0 such that 
|| M (Xl) — u i.X2)\\ H l (D) < ||Xl — X2 || £!(£>)) Vxi, X2 £ S(D), 


□ 

(35) 


where 


S(D) := {xn : Cl C D is measurable and xn(l — Xn) = 0 a.e. in D} . 

With this result it is easy to see that t m* is in fact continuous. 

Lemma 4. Let Assumptions 7|7 to 111 hold. Then the mapping t i— >■ vf := 'Ft(tq) is 
continuous from the right in 0, i.e., 

= 0 . 


lim II vf — u\ 
t\o 


HpD) 


If in addition Assumption [if^] is satisfied, then there are constants c, r > 0 such that 

||vf — u||ifi(£)) < tc, for all t E [0, t\. 

Proof. Since || • \\h} s (D) and the L 2 -norm of the gradient are equivalent norms on H t \ (D) 
it suffices to show lim t \o ||Vm* — R d) = 0. First of all introduce 

A t (x) 

which satisfies \A t {x)~ l \ < c for all t E [0, r] and hence 


Vx Ett,VtE [0, r], VC G R 2 : |C| 2 < KO^T 1 ! \£(t)(D$ t (x))- T { • (D^(x))~ T (\ 2 

= cA t (x)C • C- 


(36) 


Therefore for all f 0 E H l (D) and all t E [0,r] 

[ |V(/ 0 o 7} _1 )| 2 dx = [ A t Vf 0 -Vf 0 dx>- [ |V/ 0 | 2 dx. 

J D J D c J D 

Further, we get from this estimate that for all t E [0, r] 

c ll/o||iTi(D) < \\fo°T t ^lirqD)- (37) 

Now setting Xi := Xci and X2 := Xsu = Xn 0 Ttl ' and denoting the corresponding solutions 
of ^ by u := u(xn) and u t := w(xn t ), we infer from (|35j) and ([37]) 

c||V(u o Tt) — Vw*|| L 2 (D,n d ) < \\^ u t ~ V u \\L 2 (D,n d ) < C’llxn — Xn ° T t 1 ||l 1 (d)> 
where v* := u t o T t . Employing the previous estimates, we get for all t E [0, r] 

II Vw — VtC|| L 2 (D,R d ) < || Vu - V(«o T t ) ||i2p )R d) + ||V(m o Tt) — V?/|| L 2 (D,n d ) 

< c (||Vw — V(m o Tt)\\ L 2 ( DjR d) + ||xn ~ Xn°T t 1 ||li(d)) • 

Finally, taking into account Lemma [3| we obtain the desired continuity. The Lipschitz 
continuity under the additional Hypothesis [Ip| was shown in [JT] • □ 
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Using the previous lemma we are able to show the following. 


Lemma 5. The solution p l of (32) converges weakly in Hq(D) to the solution p of the 
adjoint equation (23)-(24). 

Proof. The existence of a solution of (32) follows from Theorem [4j Inserting if = // as 
test function in (32), we see that the estimate ||« t ||.ffi(.D) < C implies ||p*||.ffi(D) < C for 
t sufficient small. From the boundedness, we infer that (p*) t > 0 converges weakly to some 
w G Hq(D). In Lemma [4] we proved uf —> u in Hq(D ) which we can use to pass to the 
limit in rt32|) and obtain 


P 


pk 


p in Hq(D) for tk —» 0 as k —>• oo, 


where p G Hq(D) solves the adjoint equation (23)-(24). By uniqueness we conclude 
w = p. □ 


Now we proceed to the differentiation of (30) at t > 0. Introduce the notations 

:= £(t)M(t) T M(t), we obtain 


P(t) = and 

d t <5(t,p,if) = [ C(t)\B r (p) - B d oT t \ 2 ds 
J r 0 

- 2 f £, T (t)(B r (ip) - B d o T t )X7B d o T t ■ V ds 

Jr 0 

+ 2 [ d c p n (x, |M(t)V(p| 2 )(M , (t)V(p • M(t)V^)M(t)V^ • M(t)V^(t) dx 

Jd 


(38) 


Id 


/3q(x, |M(t)V^| )Q'(t)V(p ■ Vij dx 


- / /o ° T t ip£'(t) + V/o o T t • U ip£(t) dx 


id 


/ F\ty Vxp • Mi dx P f _P'(t) • M 2 dx 

^ ^ma ^^\^ma 


and this shows that for fixed ip G Hq(D ) the mapping (t, H> ip, is weakly 

continuous. This finishes proving that (H3) is satisfied. 

Consequently, we may apply Theorem[3]and obtain dJ(M V ) = 6^0(0. u,p), where u G 
Hq (D) solves the state equation 0-0 and p G H} } (D) is a solution of the adjoint equation 
(p3])-(p4l). In order to compute 9 t 0(O ,u,p), note that the integrals on T 0 vanish due to 
V = 0 on T 0 , M'(0) = - DV T , P'(0) = (div V)I 2 - DV T , Q'(O) = (div V)I 2 - DV T - DV. 

Finally we have proved formula (|27|). This concludes the proof of Theorem [5j ■ 


6 Optimization of the rotor core 

In this section, we use the shape derivative derived in Theorem[5]to determine the optimal 
design for the electric motor described in Section [2| Recall that the problem consists in 
finding the shape 0 C 11/ of the ferromagnetic subdomain of the electric motor depicted 
in Fig. [I] which minimizes the cost functional 

J(U) = [ \B r (u Q ) - B d \ 2 ds 
J r 0 
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among all admissible shapes G O where T 0 is a circle in the air gap, B t {uq) denotes 
the radial component of the magnetic flux density B = B(uq) on F 0 and B ( j is a given 
sine curve, Bd(9) = |sin(40) where 6 denotes the angle in polar coordinates with origin 
at the center of the motor; see Fig. [l] Minimizing this functional leads to a reduction of 
the total harmonic distortion (THD; see bed of the flux density which causes the rotor 
to rotate more smoothly. Here, u q is the solution of the two-dimensional magnetostatic 
boundary value problem the weak form of which reads as follows: 

Find u G Hq(D) such that / i^(|Vu|)Vw • Vndx = (/, v) for all v G Hq(D). (39) 

J D 

Here, the right-hand side / corresponds to the weak form of 

/ = /o + divM with M = MiXn m S x ) + M 2 XD\n ma ( x ) 

as in Section |H| where / 0 = J t , Mi = and M 2 = 0. The vector denotes the 

permanent magnetization of the magnets. It is a constant vector in each of the magnet 
subdomains pointing in the directions indicated in Fig. [I] and vanishes outside the magnet 
areas. Denoting M 1 - = the right hand side / G H l (D) reads 

(f,v) = f JiV + M 1 - • Vndx. (40) 

Jd 

The function Jj represents the impressed currents in the coil areas (light blue areas in 
Fig.0 and is assumed to vanish in this special application, i.e., J* = 0. 

We consider admissible shapes Q G O as in (|5]) and in Section [5} Furthermore, 

v(\ V«|) = Xn f H\ V«|) + Xd\TTJ "o 

denotes the magnetic reluctivity composed of a nonlinear function z> depending on the 
magnitude of the magnetic flux density \B\ = Vu inside the ferromagnetic material and 
a constant /y = 10 7 /(4 tt), which is expressed in the unit mkg A~ 2 s ~ 2 , otherwise. The 
constant vq is the magnetic reluctivity of vacuum which is practically the same as that 
of air. The nonlinear function v is in practice obtained from measurements and is not 
available in a closed form. However, the physical properties of magnetic fields yield the 
following characteristics of z>: 

(i) v is continuously differentiable on (0, oo), 

(ii) 3m > 0 : z>(s) > m for all s G Rq , 

(iii) P(s) < u 0 for all s G Rq , 

(iv) (P(s)s)' = P(s) + iy(s)s > m > 0, 

(v) s i —y P(s)s is strongly monotone with monotonicity constant m, i.e., 

(. P(s ) s — P(t) t ) (s — t) > m(s — t ) 2 for all s, t > 0, 

(vi) s 1 —y P(s)s is Lipschitz continuous with Lipschitz constant zy,, i.e., 

|z>(s) s — P(t) f)| < z^ol 5 — f° r all Vs,t > 0. 
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Figure 2: (a) Magnetic reluctivity 0 as a function of the magnitude \B\ = |V?/| of the 
magnetic flux density, (b) Ferromagnetic material Qf with design subdomains Q C Qf 
(highlighted). 


For more details on properties and practical realization of the function z> from given 
measurement data, we refer the reader to [16, 2"4, 35]. 

In order to be able to apply Theorem[5]to the problem above, we have to check whether 
Assumption [l] is satisfied for 

/MO := Hy/C) and MO ■= An 

Clearly, all four conditions of Assumption [I] are fulfilled for /%(() = zy, = const. Now 
let us investigate more closely M Notice the relations M|p| 2 ) = MpI) and MH 2 ) = 

MH)/(2|pI). 


1. As mentioned above, the function v is bounded from above by the magnetic reluc¬ 
tivity of vacuum uq and from below by a positive constant m, compare Fig. [sf'a). 

2. Assumption TJ2 holds for the function v by virtue of properties (v) and (vi). 

3. The numerical realization of the function 0 consists in an interproximation of given 
measurement data. The interproximation was done using splines of class C 1 . 

4. It is easy to see that this condition for /A is equivalent to 

3A, A > 0 : AIM < r] T (M|p| 2 )/ 2 + 2 MH 2 )p ®p)v< A|??| 2 , 

or in terms of z>, 


3A, A > 0 : AM 2 < V T [u(\p\)I 2 + 


\p\ 


p® p)rj< A\rj\ 2 , 
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where / 2 denotes the identity matrix of dimension 2. The eigenvalues and corre¬ 
sponding eigenvectors of the 2x2 matrix h(|p|)/ 2 + - 0 p are given by 


Ai=z>(|p|) 61 = ( (h) ’ 

a 2 = H\p\) + ^(IpDIpI v 2 = . 


From the physical properties (ii) and (iv) of z> it follows that both Ai and A 2 are pos¬ 
itive. Therefore, the assumption holds with A = min{Ai, A 2 } and A = max{Ai, A 2 }. 


Properties (v) and (vi) together with Assumption [Ip] yield existence and uniqueness of 
a solution u G Hi (D) to the state equation. Assumption [T|4| ensures the existence of an 
adjoint state p G Hq(D). 

Thus we can apply Theorem [ 5 ] and the shape derivative is given by (27). 


6.1 Numerical method 

In each iteration of the optimization process we use the shape derivative dJ{i 2; V) derived 
in ( |27| ) to compute a vector field V that ensures a decrease of the objective function 
dJ(Q]V) < 0 by displacing the interface between the iron subdomain Q and the air 
subdomain Q ref \ Q along that vector field. 


6.1.1 Setup of interface 

Due to practical restrictions we choose not to move the interface by moving the single 
points of the mesh, as it is common practice in shape optimization. Instead we model the 
interface by setting up a polygon with 151 points around each of the design subdomains 
(see Fig. [ 3 ]) and move the points of these polygons along the calculated velocity field 
V in the course of the optimization process. Each element of the design area whose center 
of gravity is inside this polygon is considered to be ferromagnetic material, the others are 
considered to be air. That way, we can avoid problems such as deformation of the fixed 
parts of the motor, i.e. magnets or the air gap, or self-intersections of the mesh. 


6.1.2 Descent Direction 

In order to get a descent in the cost functional, we compute the velocity field as follows. 
We choose a symmetric and positive definite bilinear form 


b : Hf\ (D rot ) x A/]] ( D rot ) —>• M 


defined on the subdomain D rot of D representing the rotor and compute V as the solution 
of the variational problem: 

Find V E P h such that b(V, W) = -dJ(Q, W) for all IV G P h , (41) 

where P/, C H ( \ (D rot ) is a finite dimensional subspace. Outside D rot we extend V by 
zero. Note that, by this choice, the condition V = 0 on r 0 , which is assumed in Section 
[5| is satisfied. The obtained descent directions V G Ph will also be in W 1,oc (D, R 2 ) 
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a) 

Figure 3: (a) Interface points around eight design areas of the electric motor, (b) Zoom 
on the upper design area: Fictitious interface polygon consists of 151 points (71 on the 
left, 71 on the right and 9 on top) 


and, consequently, they are admissible vector fields defining a flow 7)' . The solution V 
computed this way is a descent direction for the cost functional since 


dJ(Q, V ) = -b(V, V) < 0. 


For our numerical experiments, we choose the bilinear form 
b : Hq (D mt ) x H Xj (D roX ) —> M 

b(V,W)= [ (a DV : DW + V ■ W) dx. 


'D r, 


Here, the penalization function a E L°°(D rot ) is chosen as 

1 x E Q 

a(x) = ^ 10 x E fl £ \fl 
10 2 else, 


(42) 


where fl £ = {x E D rot : distfa;, fl) < z\ for some small e > 0. With this choice of a , we 
ensure that the resulting velocity field V is small outside the design region O ref . 

For all numerical simulations, we used piecewise linear finite elements on a triangular 
grid with 44810 degrees of freedom and 89454 elements where we chose a particularly 
fine discretization in the design regions Q rcf (53488 design elements). The nonlinear state 
equation © is solved by Newton’s method. All arising linear systems of finite element 
equations are solved using the parallel direct solver PARDISO 


6.1.3 Updating the interface 

For updating the interface, we perform a backtracking (line search) algorithm: Once a 
descent direction V is computed, we move all interface points a step size r = T vnit in the 
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direction given by V and evaluate the cost function for the updated geometry. If the cost 
value has not decreased, the step size r is halved and the cost function is evaluated for 
the new, updated geometry. We repeat this step until a decrease of the cost function has 
been achieved. When the step size becomes too small such that no element switches its 
state, the algorithm is stopped. 


6.2 Numerical Results 

The procedure is summarized in Algorithm [lj 

Algorithm 1. Set converged = false 
While !converged 

1. Compute state u using Q and adjoint state p using (22). 


2. Compute shape derivative dJ(kl;V) from (27). 


3. Compute descent direction V using (41). 

4■ Find step width r that yields a decrease in the cost function using backtracking. 
5. If decrease in the cost function could be achieved: 

Update interface and go to step 1. 

else: 


converged = true. 

The final design after 35 iterations of Algorithm [l] can be seen in Fig. [4} The cost 
function is reduced from 1.3033 * 10 -3 to 0.94997 * 10 —3 , i.e., by about 27%. The radial 
component of the magnetic field on the circle T 0 in the air gap for the initial (blue), 
desired (green) and final design (red) can be seen in Fig. [5j The optimization process 
took about 26 minutes using a single core on a laptop. 


7 Conclusion 

In this paper, we have performed the rigorous analysis of the shape sensitivity analysis 
of a subregion Q of the rotor of an electric motor in order to match a certain rotation 
pattern. The shape derivative of the cost functional was computed efficiently using a 
shape-Lagrangian method for nonlinear partial differential equation constraints which 
allows to bypass the computation of the material derivative of the state. The implemen¬ 
tation of the obtained shape derivative in a numerical algorithm provides an interesting 
shape which allows us to improve the rotation pattern. In the numerical experiment 
presented in this work, we chose a rather simple way of updating the interface: We just 
switched the state of single elements of the finite element mesh as discrete entities. For a 
more accurate resolution of the interface, one may employ a nonstandard finite element 
method like the extended Finite Element Method (XFEM) [33j or the immersed FEM 
[30j, or a discontinuous Galerkin approach based on Nitsche’s idea, see [13]. These ap¬ 
proaches make it possible to represent an interface that is not aligned with the underlying 


21 






Figure 4: (a) Final design of one component 0/,, after 35 iterations of Algorithm [TJ (b) 
Upper-right quarter of the optimized motor. 



Figure 5: (a) Radial component of the magnetic flux density B along curve r 0 in air gap: 
initial design (blue), desired sine curve B^ (green), final design (red). (b) Decrease of cost 
functional in the course of the optimization process. 
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FE discretization without loss of accuracy. An alternative way to achieve this would be 
to locally modify the finite element basis in a parametric way as it is done in (8j. 
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